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Reduced-order modeling (ROM) methods are applied to the CFD-based aeroelastic anal- 
ysis of the AGARD 445.6 wing in order to gain insight regarding well-known discrepancies 
between the aeroelastic analyses and the experimental results. The results presented in- 
clude aeroelastic solutions using the inviscid CAP-TSD code and the FUN3D code (Euler 
and Navier-Stokes). Full CFD aeroelastic solutions and ROM aeroelastic solutions, com- 
puted at several Mach numbers, are presented in the form of root locus plots in order to 
better reveal the aeroelastic root migrations with increasing dynamic pressure. Impor- 
tant conclusions are drawn from these results including the ability of the linear CAP-TSD 
code to accurately predict the entire experimental flutter boundary (repeat of analyses 
performed in the 1980’s), that the Euler solutions at supersonic conditions indicate that 
the third mode is always unstable, and that the FUN3D Navier-Stokes solutions stabilize 
the unstable third mode seen in the Euler solutions. 


I. Introduction 

Classical linear aeroelastic analyses typically produce Velocity- Damping- Frequency (V-g-f) plots and/or 
root locus plots. The use of these plots has enabled the aeroelastician to view the nature of the flutter 
mechanism(s) in addition to identifying the condition(s) at which flutter occurs. The rapid creation of these 
plots was facilitated by the use of linear unsteady aerodynamics and linear aeroelastic equations of motion . * 1 

During the last few years, higher-order CFD-based methods have become an important method for the 
computation of nonlinear unsteady aerodynamics for use in aeroelastic analyses. The use of these higher- 
order methods provides valuable insight regarding complex flow physics at conditions where linear methods 
are not theoretically valid. However, the increased computational cost associated with the computation of 
unsteady aerodynamics and aeroelastic responses using higher-order methods has resulted in a subtle change 
in the manner in which the aeroelastician evaluates and interprets these analyses. First, the increased 
computational cost of these analyses has tended to dictate a ’’snapshot” approach to aeroelastic analyses 
whereby the aeroelastic response at a handful of dynamic pressures is all that is computed. This ’’snapshot” 
approach is used to identify the flutter dynamic pressure but the actual flutter mechanism is not easily 
discernible. Second, due to the complexity of the computational methods, methods that could rapidly 
generate V-g-f plots and/or root locus plots were not available. However, with the development of reduced- 
order modeling (ROM) methods , 2-4 the rapid generation of root locus plots using CFD-based unsteady 
aerodynamics is now available to aeroelasticians. 

The goal behind the development of a ROM for the rapid computation of unsteady aerodynamic and 
aeroelastic responses is aimed at addressing two challenges. The first challenge is the computational cost as- 
sociated with full CFD aeroelastic simulations, which increases with the fidelity of the nonlinear aerodynamic 
equations to be solved as well as the complexity of the configuration. Computational cost, however, may 
be reduced via the implementation of parallel processing techniques, advanced algorithms, and improved 
computer hardware processing speeds. 

The second, more serious, challenge is that the information generated by these simulations cannot be 
used effectively within a preliminary design environment. For this reason, parametric variations and design 
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studies can only be performed by trial-and-error. As a result, the integration of computational aeroelastic 
simulations into preliminary design activities involving disciplines such as aeroelasticity, aeroservoelasticity 
(ASE), and optimization continues to be a costly and impractical venture. 

Development of a ROM entails the development of a simplified mathematical model that captures the 
dominant dynamics of the original system. This alternative mathematical representation of the original 
system is, by design, in a mathematical form suitable for use in a multidisciplinary, preliminary design envi- 
ronment. As a result, use of the ROM in other disciplines is possible, thereby addressing the second challenge. 
The simplicity of the ROM yields significant improvements in computational efficiency as compared to the 
original system, thereby, also addressing the first challenge. 

A CFD-based model of an aeroelastic system (such as the FUN3D code 5,6 ) consists of the coupling 
of a nonlinear unsteady aerodynamic system (flow solver) with a linear structural system as depicted in 
Figure 1. Standard CFD-based aeroelastic analyses are performed via iterations between the nonlinear 
unsteady aerodynamic system and the linear structural system. Throughout this paper, these standard 
CFD-based solutions will be referred to as ’’Full CFD Solutions”. 

The present study also involves the generation of linearized unsteady aerodynamic ROMs (in state- 
space form), using the general procedure depicted in Figure 2. In this situation, the linear structural system 
within the CFD code is bypassed so that only the nonlinear unsteady aerodynamic system is excited. Specific 
modal inputs are applied to the nonlinear unsteady aerodynamic system and the generalized aerodynamic 
force (GAF) outputs from this solution, along with the inputs, are used in a system identification process 
to create the linearized unsteady aerodynamic ROM. This unsteady aerodynamic ROM is then coupled to 
a state-space model of the structure in order to create the aeroelastic simulation ROM. The aeroelastic 
simulation ROM is then used for aeroelastic analyses. For the discussions that follow, the term ROM will 
refer to the unsteady aerodynamic state-space model. When the unsteady aerodynamic state-space model 
(ROM) is connected to a state- space model of the structure, this system is often also referred to as a ROM. 
However, to avoid confusion, the aeroelastic system consisting of an unsteady aerodynamic ROM (state-space 
form) coupled with a linear modal model of the structure (also state-space form) will be referred to as the 
aeroelastic simulation ROM. 



Figure 1. Coupling of linear structural model and nonlinear unsteady aerodynamics within an aeroelastic CFD code 
such as FUN3D. 
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Figure 2. Generation of generalized aerodynamic forces (GAFs) used for system identification process. 


The development of CFD-based ROMs continues to be an area of active research at several national 
and international government, industry, and academic institutions. 7-11 Development of ROMs based on 
the Volterra theory is one of several ROM methods that has received attention over the last few years. 12 19 
Although the primary focus of this paper is the development and application of unsteady aerodynamic ROMs 
for subsequent use in aeroelastic analyses, the development of ROMs for the rapid computation of nonlinear 
stability and control derivatives using CFD codes 20, 21 is an active area of research as well. 

Silva and Bartels 8 introduced the development of linearized, unsteady aerodynamic state-space models 
for prediction of flutter and aeroelastic response using the parallelized, aeroelastic capability of the CFL3Dv6 
code. The results presented provided an important validation of the various phases of the ROM development 
process. The Eigensystem Realization Algorithm (ERA), 22 which transforms an impulse response (one 
form of a ROM) into state-space form (another form of a ROM), was applied to the development of the 
aerodynamic state-space models. The ERA is part of the SOCIT (System/Observer/Controller Identification 
Toolbox). 23 Flutter results for the AGARD 445.6 aeroelastic wing using the CFL3Dv6 code were presented, 
including computational costs. 8 Unsteady aerodynamic state-space models were generated and coupled 
with a structural model within a MATLAB/SIMULINK 24 environment for rapid calculation of aeroelastic 
responses including the prediction of flutter. Aeroelastic responses computed directly using the aeroelastic 
simulation ROM showed excellent comparison with the aeroelastic responses computed using the CFL3Dv6 
code. 25 

Previously, 8 the aerodynamic impulse responses (unit pulses) that were used to generate the unsteady 
aerodynamic state-space model were computed via the excitation of one mode at a time. For a four-mode 
system such as the AGARD 445.6 wing, these computations are not very expensive. However, for more 
realistic cases where the number of modes can be an order of magnitude or more larger, the one- mode- at-a- 
time method becomes prohibitively expensive. Towards the solution of this problem, new methods have been 
developed. Kim et al. 10 have proposed methods that enable the simultaneous application of structural modes 
as CFD input, greatly reducing the cost of identifying the aerodynamic impulse responses from the CFD 
code. Kim’s method consists of using simultaneous staggered step inputs, one per mode, and then recovering 
the individual responses from this simultaneous excitation. Silva 2 has developed a method that enables the 
simultaneous excitation of the structural modes using orthogonal functions. Both of these methods require 
only a single CFD solution and the methods are independent of the number of structural modes. Silva 26 has 
also developed a method for generating static aeroelastic solutions and matched-point aeroelastic solutions 
using a ROM. The methods developed by Silva 2,26 have already been implemented in the FUN3D CFD 
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code. In addition, methods for generating root locus plots of the combined structural state-space model 
and unsteady aerodynamic state-space model were developed by Silva et al. 3 These ROM-based root locus 
methods were applied to fixed-wing configurations and subsequently to launch vehicle configurations. 4 The 
present paper will focus on the application of these ROM and root locus methods in order to visualize the 
aeroelastic behavior of the AGARD 445.6 wing as a function of Mach number, dynamic pressure, and fluid 
dynamic equation (CAP-TSD, inviscid FUN3D, and viscous FUN3D). 

The paper begins with a description of the AGARD 445.6 wing and a comparison of experimental and 
computational flutter results obtained to date by various researchers. 27 Computational methods and related 
models are introduced including the CAP-TSD code, the FUN3D code (inviscid and viscous grids), and 
the FUN3D ROM creation process. The results to be presented are grouped into two categories. The first 
category consists of the ’’Full CFD Solutions” based on the standard iterative approach briefly described 
above. These results will include full CAP-TSD solutions and full FUN3D solutions (inviscid and viscous) 
for several Mach numbers and dynamic pressures. The aeroelastic transients computed via the full CFD 
solutions are analyzed for their damping and frequency content in order to generate aeroelastic root locus 
plots. The second category of results consists of the ”FUN3D ROM Solutions”. At present, the ROM 
method has not been implemented in the CAP-TSD code, therefore all ROM solutions will be FUN3D ROM 
solutions. The FUN3D ROM solutions will be presented in the form of root locus plots at several Mach 
numbers. Finally, some concluding remarks will be provided. 

II. AGARD 445.6 Wing 

The AGARD 445.6 wing was tested in the NASA Langley Transonic Dynamics Tunnel (TDT) in 1961. 28 
Flutter data from this test has been publicly available for over 20 years and has been widely used for 
preliminary computational aeroelastic benchmarking. The AGARD wing planform was sidewall-mounted 
and had a quarter-chord sweep angle of 45 degrees, an aspect ratio of 1.65, a taper ratio of 0.66, a wing 
semispan of 2.5 feet, a wing root chord of 1.833 feet, and a symmetric airfoil. The wing was flutter tested 
in both air and R-12 heavy gas test mediums at Mach numbers from 0.34 to 1.14 at zero-degrees angle 
of attack. Unfortunately, this data set lacks unsteady surface pressure measurements necessary for more 
extensive code validation. 

A broad range of FUN3D computations 27 for the AGARD 445.6 wing were performed across the entire 
Mach number range of the experimental data, assuming both inviscid and viscous flows with air as the 
working fluid. The first four structural modes were used in the aeroelastic analysis and are presented in 
Figure 3. 

Figure 4 and Figure 5 present comparisons among the experimental flutter speed index and frequency ratio 
values, respectively, with those obtained using the FUN3D code, and those published in the literature. 29-32 
In general, in the subsonic flow regime, the computational data matches the experimental data well, while 
a broad range in the computational data is observed in the high subsonic and supersonic flow regimes. 
Chwalowski et al 27 emphasize the importance of applying the viscous flow assumption (using FUN3D) at 
the high subsonic conditions as evidenced by the improved result over the FUN3D inviscid subsonic solutions. 
In addition, Chwalowski et al 27 carried out a grid refinement study that improved the correlation between the 
viscous FUN3D solution and the experiment at the supersonic Mach numbers. However, an important point 
to be made is that the FUN3D (inviscid and viscous) flutter solutions at the supersonic Mach numbers are 
for a flutter mechanism comprised of the coalescence of the first and second modes. Based on the frequency 
ratio of the results presented from other references, it appears that those results correspond to the same 
mechanism. This assumption needs to be confirmed with each individual researcher that has provided the 
results for this wing included in these figures. 
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Mode 1 



(a) Mode 1 


Mode 2 



(b) Mode 2 


Mode 3 



(c) Mode 3 



(d) Mode 4 


Figure 3. The first four modes of the AGARD 445.6 wing where ’zmd’ is the modal deflection in the z-direction. 
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Figure 4. Experimental and computational flutter speed index versus Mach number for the AGARD 445.6 wing. 
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Figure 5. Experimental and computational frequency ratio versus Mach number for the AGARD 445.6 wing. 


In order to simplify the discussion regarding comparisons of experimental and various full CFD and ROM 
solutions, our results will be presented in terms of flutter dynamic pressure (psf) and flutter frequency (Hz) 
for a subset of FUN3D results only. Figure 6 presents flutter boundaries in terms of flutter dynamic pres- 
sure in psf for Experiment, FUN3D/Euler (FUN3D-E), FUN3D/Navier-Stokes/Spallart-Allmaras turbulence 
model for the baseline grid (FUN3D NS SA Baseline Grid), and FUN3D/Navier-Stokes/Spallart-Allmaras 
turbulence model for the fine grid (FUN3D NS SA Fine Grid). The comparison between the baseline grid 
and the fine grid is the result of the grid refinement study previously mentioned. Figure 7 presents the 
corresponding flutter frequencies (Hz) as a function of Mach number. All subsequent FUN3D results in 
this paper will include Euler solutions and Navier- Stokes solutions (with the SA turbulence model) for the 
baseline grid. 

As will be discussed in subsequent sections of this paper, the discrepancy between the computational 
and experimental results at the supersonic Mach numbers raises interesting questions. The first question 
concerns results obtained by Yates 33 using modified strip analysis as well as results obtained by Bennett et 
al. 34 in the late 1980’s using the CAP-TSD code (linear and non-linear solutions). Both references indicate 
that the entire flutter boundary is well predicted using linear methods. These results are consistent with 
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Figure 6. Experimental and computational flutter dynamic pressure (psf) versus Mach number for the AGARD 445.6 
wing. 


the fact that the AGARD 445.6 wing has a thin airfoil and the wing does not reach transonic conditions 
until about M=0.98. At the supersonic Mach numbers of interest (M=1.07 and M=1.14), the flow is entirely 
supersonic and should therefore be well predicted with the use of linear unsteady aerodynamic methods. 
Given these facts, the first question is: Why are the Euler and Navier- Stokes solutions so different from 
the linear solutions at conditions where there is, essentially, no nonlinearity in the flow? Furthermore, the 
flutter results computed using the linear equations within the CAP-TSD code serve to correct statements 
often found in other references that classify the flutter boundary of the AGARD 445.6 wing as a nonlinear 
transonic flutter dip. Clearly, if linear CAP-TSD can accurately predict the flutter boundary at high subsonic 
Mach numbers and if transonic flow is not present at Mach numbers below M=0.98, the flutter boundary 
up to that Mach number cannot be referred to as a transonic flutter dip. The dip in the flutter boundary 
that is observed is due, primarily, to compressibility and is not the result of a nonlinear transonic effect. 

The second question raised by the discrepancy in the computational and experimental flutter results 
at the supersonic Mach numbers is: What is the actual flutter mechanism at supersonic conditions? The 
experimental value of flutter frequency at M=1.141 is about 17.5 Hz, suggesting that the measured flutter 
mechanism involved the first and second modes. Based on this indication, it appears that computational 
results presented to date have focused on this particular flutter mechanism. However, it will be shown later 
in the paper that the third mode is always unstable for the inviscid (Euler) solutions. It is not clear if this 
situation is true for all other inviscid solutions performed by other researchers but it is certainly true for the 
inviscid solutions computed using FUN3D. The authors have confirmed that this is the case as well for the 
inviscid solution from the CFL3Dv6 code (a NASA Langley-developed structured CFD aeroelastic code), 
although those results are not presented here. The use of FUN3D ROMs at this condition served to highlight 
this issue with the obvious clarity of a root locus plot to be presented later in the paper. Although this 
discussion may raise more questions than it answers, the authors feel this is an important discussion to be 
raised if we are to seriously address the validation of our computational aeroelastic tools. 

III. Computational Methods 


III. A. CAP-TSD Code 

The CAP-TSD (Computational Aeroelasticity Program- Transonic Small Disturbance) code is a finite dif- 
ference program that solves the general- frequency modified TSD potential equation. The TSD potential 
equation is solved within CAP-TSD by a time-accurate approximate factorization (AF) algorithm developed 
by Batina. 35 The CAP-TSD program can be used for the analysis of configurations with combinations 
of lifting surfaces and bodies including canard, wing, tail, control surfaces, tip launchers, pylons, fuselage, 
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Figure 7. 


Experimental and computational flutter frequency (Hz) versus Mach number for the AGARD 445.6 wing. 


stores, and nacelles. The CAP-TSD code was applied to several configurations 34,36,37 in the late 1980’s and 
the early 1990’s. 

In the present effort, linear CAP-TSD solutions were generated, repeating the work performed in the late 
1980’s. A new contribution to that original effort was the inclusion of root locus plots generated from the 
CAP-TSD aeroelastic transient responses. The procedure for generating the aeroelastic transients at various 
Mach numbers and dynamic pressures using CAP-TSD is the same as the procedure for FUN3D described 
in a subsequent section. 

III.B. FUN3D Code, Grids, and Analysis Procedure 

The following subsections describe the parallelized, aeroelastic version of the unstructured mesh solver 
FUN3D code, the inviscid and viscous grids used, and a brief review of what will be referred to as the 
FUN3D full solution (in contrast to the FUN3D ROM solution). 

III.B. 1. FUN3D Code 

The unstructured mesh Euler (inviscid)/Navier-Stokes (viscous) solver used for this study is FUN3D. 5 Within 
the code, the unsteady Navier-Stokes equations are discretized over the median dual volume surrounding each 
mesh point, balancing the time rate of change of the averaged conserved variables in each dual volume with 
the flux of mass, momentum and energy through the instantaneous surface of the control volume. Additional 
details regarding the aeroelastic capability within the FUN3D code can be found in the references. 38,39 

III.B. 2. FUN3D Grids 

Unstructured tetrahedral grids used in this study were generated using VGRID 40 with input prepared using 
GridTool. 41 For the AGARD 445.6 wing grids, the boundary layer consisted of tetrahedral elements. Only 
two grids were used for the present analyses: an inviscid (Euler) grid consisting of two million nodes and a 
viscous (Navier-Stokes) grid consisting of four million nodes also referred to as the baseline grid. Figure 8 
shows the planform, surface grids, and the surface C p at M= 1.141 for the inviscid and viscous grids used 
in this analysis. A relatively strong compression near the trailing edge is present in the inviscid solution. 
Although the flow is entirely supersonic, it is not surprising to see a strong compression as a result of the 
inviscid (Euler) analysis. 

Note that the surface grid for both the inviscid and viscous grids is the same with the difference in grid 
dimensions accounting for a denser grid normal to the surface for the viscous grid (not visible in this figure). 
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The surface C p at M= 1.141 for the viscous solution indicates that the inclusion of viscosity has diminished 
the strong compression seen for the inviscid solution. 




(a) Inviscid FUN3D surface grid for the 
445.6 wing. 


(b) Pressure coefficient contours on the 
surface for the inviscid FUN3D solution. 




(c) Viscous surface grid for the 445.6 
wing. 


(d) Pressure coefficient contours on the 
surface for the viscous FUN3D solution. 


Figure 8. Inviscid and viscous grids and pressure distributions of the AGARD 445.6 wing at M=1.141. 


III.B.3. FUN3D Analysis Procedure 

Solutions to the Reynolds- Aver aged Navier-Stokes (RANS) equations were computed using the FUN3D flow 
solver. The FUN3D solutions presented in this paper were obtained with an augmented van Leer limiter, 42 
Low-Difussion Flux-Splitting Scheme (LDFSS) 43 for inviscid fluxes, and the Spalart-Allmaras turbulence 
model. 44 For the asymptotically steady cases under consideration, time integration was accomplished by an 
Euler implicit backwards difference scheme, with local time stepping to accelerate convergence. Most of the 
steady state cases in this study were run for 5000 iterations to achieve convergence of forces and moments 
to within ±0.5% of the average of their last 1000 iterations. 

In order to perform static and dynamic aeroelastic solutions, interpolation of the structural mode shapes 
onto the CFD surface grid is required. This interpolation is done as a preprocessing step. 45 The final surface 
deformation at each time step is a linear superposition of all the modal deflections. 

The standard procedure for performing a FUN3D full aeroelastic solution was performed as follows. First, 
the steady CFD solution was obtained on the rigid vehicle. Next, a static aeroelastic solution was obtained 
by continuing the CFD analysis in a time accurate mode, allowing the structure to deform but with a high 
value of structural damping (0.99) so the structure could find its equilibrium position with respect to the 
mean flow before the dynamic response was started. However, for the AGARD 445.6 wing, there was no 
need to compute static aeroelastic results due to the fact that the airfoil is symmetric and the analyses are 
performed with the wing at zero degrees angle of attack. Finally, for the dynamic response, the damping 
was set to an assumed value (0.00), and the structure was perturbed in generalized velocity for each of the 
four modes included in the structural model. The flow was then solved in the time accurate mode. 
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III.C. FUN3D Reduced Order Models 


In this section, the FUN3D reduced order model (ROM) development process is briefly reviewed with the 
details being deferred to the references. 

An outline of the FUN3D ROM development process is as follows: 

1. Generate the number of functions from a selected family of orthogonal functions 2 that corresponds to 
the number of structural modes; 

2. Apply the generated input functions simultaneously via one CFD execution resulting in GAF responses 
due to these inputs; these responses are computed directly from the restart of a steady rigid CFD solution 
(not about a particular dynamic pressure); 

3. Using the simultaneous input/output responses, identify the individual impulse responses using the 
PULSE algorithm within SOCIT; 46 

4. Transform the individual impulse responses generated in Step 3 into an unsteady aerodynamic state- 
space system using the ERA 22 (within SOCIT); 

5. Evaluate/ validate the state-space models generated in Step 4 via comparison with CFD results (i.e. , 
ROM results vs. full CFD solution results); 

A schematic of steps 1-4 of the process outlined above is presented as Figure 9. Using modal information 
(generalized mass, frequencies, and dampings), a state-space model of the structure is generated. This 
state-space model of the structure is referred to as the structural state-space ROM (Figure 10). Once an 
unsteady aerodynamic ROM and a structural state-space ROM have been generated, they are combined 
to form an aeroelastic simulation ROM (see Figure 11). Then root locus plots are extracted from the 
aeroelastic simulation ROM. 
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Figure 9. Improved process for generation of an unsteady aerodynamic ROM (Steps 1-4). 
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Figure 10. Process for generation of a structural state-space ROM. 
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Figure 11. Process for generation of an aeroelastic simulation ROM consisting of an unsteady aerodynamic ROM and 
a structural state-space ROM. 


IV. Results 


IV. A. Linear CAP-TSD Results 

In this section, linear CAP-TSD results are presented at three Mach numbers: M=0.90, 0.96, and 1.141. 
Results are presented in the form of aeroelastic transients for different dynamic pressures and a root locus 
plot per Mach number. The ROM method has not been incorporated into the CAP-TSD code at this point 
in time. Instead, the aeroelastic transients are analyzed for damping and frequency content using a new 
MATLAB-based version of a procedure originally developed by Bennett and Desmarais. 47 This algorithm 
applies MATLAB’s curve fitting routines along with MATLAB’s optimization routines in order to find the 
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best curve fit for a given aeroelastic transient consisting of amplitude, damping, and frequency for up to 
three combined sinusoidal functions. 

For the sake of brevity, a small number of aeroelastic transients are analyzed per Mach number using this 
new algorithm. For example, four dynamic pressures were analyzed for damping and frequency content for the 
M=0.90 results, three dynamic pressures were analyzed for damping and frequency content for the M=0.96 
results, and five dynamic pressures were analyzed for damping and frequency content for the M=1.141 results. 
Additional damping and frequency estimates could be generated at additional dynamic pressures per Mach 
number in order to generate denser root locus plots. However, as long as the primary flutter mechanism was 
evident, a need for additional damping and frequency estimates was deemed unnecessary. 

IV.A.l. Mach number 0.9 

Figure 12 presents the linear CAP-TSD generalized coordinates for all four modes at M=0.90 and a dynamic 
pressure of 50 psf. As is clear from the figure, this condition is stable. 



Figure 12. Linear CAP-TSD generalized coordinates at M=0.90, Q=50 psf 


Figure 13 presents the linear CAP-TSD generalized coordinates for all four modes at M=0.90 and a 
dynamic pressure of 100 psf. At this condition an unstable flutter condition is evident. Figure 14 presents 
the aeroelastic root locus as a function of dynamic pressure for M=0.90. This root locus plot indicates that 
flutter occurs just above 90 psf and that the flutter mechanism is dominated by the first mode participation. 
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Figure 13. Linear CAP-TSD generalized coordinates at M=0.90, Q=100 psf 


IV. A. 2. Mach number 0.96 

Figure 15 presents the linear CAP-TSD generalized coordinates for all four modes at M=0.96 and a dynamic 
pressure of 50 psf. As is clear from the figure, this condition is stable. Figure 16 presents the linear CAP-TSD 
generalized coordinates for all four modes at M=0.96 and a dynamic pressure of 90 psf. At this condition an 
unstable flutter condition is evident. Figure 17 presents the aeroelastic root locus as a function of dynamic 
pressure for M=0.96. This root locus plot indicates that flutter occurs just above 75 psf and that the 
flutter mechanism is, again, dominated by the first mode participation. For this case, the MATLAB-based 
damped sine curve-fitting function had some difficulty in providing good estimates of damping for an unstable 
transient. 
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Figure 14. Linear CAP-TSD aeroelastic root locus for M=0.90. 



Figure 15. Linear CAP-TSD generalized coordinates at M=0.96, Q=50 psf 
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Figure 16. Linear CAP-TSD generalized coordinates at M=0.96, Q=90 psf 



Figure 17. Linear CAP-TSD aeroelastic root locus for M=0.96 
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IV. A. 3. Mach number 1.1 1^1 


Figure 18 presents the linear CAP-TSD generalized coordinates for all four modes at M=1.141 and a dynamic 
pressure of 50 psf. As is clear from the figure, this condition is stable. Figure 19 presents the linear CAP-TSD 
generalized coordinates for all four modes at M=1.141 and a dynamic pressure of 100 psf. This condition 
is stable as well. Figure 20 presents the linear CAP-TSD generalized coordinates for all four modes at 
M=1.141 and a dynamic pressure of 140 psf. Clearly, at this condition, an unstable flutter condition is 
evident. Figure 21 presents the aeroelastic root locus as a function of dynamic pressure for M= 1.141. This 
root locus plot indicates that flutter occurs at about 140 psf. Once again, the flutter mechanism is dominated 
by the first mode participation. 



Figure 18. Linear CAP-TSD generalized coordinates at M=1.141, Q=50 psf 


16 of 29 


American Institute of Aeronautics and Astronautics 



Figure 19. Linear CAP-TSD generalized coordinates at M=1.141, Q=100 psf 



Figure 20. Linear CAP-TSD generalized coordinates at M=1.141, Q=140 psf 
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Figure 21. Linear CAP-TSD aeroelastic root locus for M=1.141 


Flutter dynamic pressures previously shown in Figure 6, but now modified to include the linear CAP- 
TSD results, are presented in Figure 22. As mentioned before and as presented by Cunningham et al. 36 
and Bennett et ah, 34 the linear CAP-TSD code shows excellent agreement with experiment for the subsonic 
Mach numbers. For the supersonic Mach number, the linear CAP-TSD result is above (not conservative) 
the experimental flutter dynamic pressure, as are all the FUN3D results, but it is closer to the experimental 
value than any of the FUN3D results presented. 
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Figure 22. Flutter boundaries as a function of dynamic pressure (psf) including linear CAP-TSD results. 


Flutter frequencies corresponding to the flutter dynamic pressures presented in Figure 22 are presented 
in Figure 23. For the subsonic Mach numbers, the linear CAP-TSD flutter frequencies are higher than the 
experimental and the FUN3D flutter frequencies. However, for the supersonic condition, the linear CAP-TSD 
flutter frequency is closer to the experimental value than the FUN3D results. 
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Figure 23. Flutter frequencies (Hz) including linear CAP-TSD results. 


IV. B. FUN3D Full and ROM Results 

In this section aeroelastic transients and aeroelastic root locus plots are presented for the FUN3D full and 
ROM results for both inviscid and viscous solutions. The FUN3D full results consist of aeroelastic transients 
at various dynamic pressures for two Mach numbers: M=0.96 and M= 1.141. The FUN3D ROM results will 
consist of aeroelastic root locus plots for the same Mach numbers. The root locus plots generated using 
the ROM procedure described above are used to identify aeroelastic behavior and flutter mechanisms. The 
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aeroelastic transients generated using the FUN3D full solutions are used to validate the FUN3D ROM results 
at specific dynamic pressures. 

IV.B.l. Inviscid Results 

In this section, inviscid FUN3D results are presented for both full and ROM solutions. Presented in Figure 24 
is the aeroelastic root locus plot for M=0.96 generated using the FUN3D ROM method. In contrast to the 
root locus plots presented for the linear CAP-TSD solutions, these root locus plots contain root values at 
twenty dynamic pressures from zero to 114 psf. The reason for this increased resolution in dynamic pressure 
values is the ROM procedure itself. In this case, there is no need to analyze an aeroelastic transient at each 
and every dynamic pressure as was the case for the linear CAP-TSD results. Because the ROM procedure 
generates a combined aeroelastic state-space model that consists of a state-space model of the structure and 
a state-space model of the unsteady aerodynamics (from FUN3D), root locus plots can be generated for 
any number and any increment of dynamic pressure in a matter of seconds. The flutter mechanism for the 
FUN3D inviscid solution at this Mach number is clearly dominated by the first mode with some coupling 
with the second mode noticeable. The third and fourth modes are stable. 
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Figure 24. FUN3D ROM aeroelastic root locus plot for M=0.96, inviscid solution. 


In order to better visualize the root migrations for the first mode, a zoomed-in version of the root locus 
plot is presented as Figure 25. The increment in dynamic pressure for this root locus plot is 6 psf starting 
with 0 psf yielding a flutter dynamic pressure of approximately 30 psf. The flutter dynamic pressure is, 
therefore, about 30 psf. This result is very close to and consistent with the FUN3D full solution flutter 
dynamic pressure presented in Figure 6. However, the inviscid result at this Mach number does not compare 
well with the experiment. This is not surprising as inviscid solutions tend to have stronger shocks that are 
farther aft and therefore induce a stronger and earlier onset of flutter. When viscosity is introduced into 
the solution, the shock strength is reduced and the shock position is moved forward resulting in the onset of 
flutter at a higher dynamic pressure. This effect is discussed in the next section of this paper. 
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Figure 25. Zoomed-in version of FUN3D ROM aeroelastic root locus plot for M=0.96, inviscid solution. 


Solutions are now presented for the supersonic Mach number of 1.141. As presented in Figure 6, there is a 
wide variation of flutter dynamic pressures and flutter frequencies at this condition. The discrepancy between 
many of the solutions and the experiment as well as the discrepancy amongst the various solutionmethods 
has long been a source of speculation. Although the authors do not present a conclusive answer to the source 
of these discrepancies, it is hoped that the results presented will spur additional discussion and research. 

Figure 26 presents the aeroelastic root locus plot for M=1.141 generated using the FUN3D ROM method. 
There are two obvious flutter mechanisms: one flutter mechanism is an instability involving the first mode 
while the other flutter mechanism involves an instability of the third mode. The third mode is always 
unstable while the first mode instability occurs at a dynamic pressure of about 300 psf. The third mode 
instability was a surprise in that it is not mentioned by other researchers. So in order to validate the accuracy 
of this aeroelastic root locus plot, the generalized coordinates from a FUN3D full solution are inspected. 

Presented in Figure 27 are the aeroelastic transients for the four modes at M= 1.141 and a dynamic 
pressure of 30 psf. The first mode, with the largest amplitude, is clearly stable. However, discerning the 
stability for the other three modes with smaller and similar amplitudes is not as obvious as it is for the first 
mode. If only the third mode is visualized, as in Figure 28, the unstable nature of this mode becomes clearer. 
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Figure 26. FUN3D ROM aeroelastic root locus plot for M=1.141, inviscid solution. 
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Figure 27. FUN3D full solution generalized coordinates at M=1.141, Q=30 psf, inviscid solution. 
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Figure 28. FUN3D full solution third generalized coordinate at M=1.141, Q=30 psf, inviscid solution. 


This third mode instability, not mentioned in any other publications on the flutter boundary of the 
AGARD 445.6 wing, raises an important question. Is this third mode instability present in all inviscid 
(Euler) solutions of the AGARD 445.6 wing presented in the literature? The answer to this question requires 
consultation with researchers that have performed inviscid CFD aeroelastic analyses for this wing at this 
condition. However, as mentioned previously, it appears that the focus of all previous inviscid analyses at 
supersonic conditions was the first mode instability. If that is the case, it is easy to see how the third mode 
instability might have been ignored. In addition, for analyses performed in the early days of computational 
aeroelasticity, Figure 27 would have consisted of fewer time steps (due to computational cost at the time), 
thereby making it difficult to visually notice the third mode instability. It should be stated that the authors 
have confirmed the existence of this third mode instability in previous solutions obtained using the CFL3D 
code. 

IV. B. 2. Viscous Results 

In this section, viscous FUN3D full and ROM solutions are presented at M= 1.141 since significant discrep- 
ancies, as just discussed, exist at this condition. The results for FUN3D full and ROM solutions at subsonic 
Mach numbers agree well with each other and with experiment and are not presented in this paper. 

Presented in Figure 29 is the root locus plot generated using the FUN3D ROM viscous solution at 
M=1.141 in dynamic pressure increments of 6 psf to 114 psf. It appears that the third mode instability 
exhibited by the inviscid solution has been stabilized by the inclusion of viscous effects. 
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Figure 29. Viscous FUN3D ROM root locus plot at M=1.141. 


Presented in Figure 30 is a zoomed- in version of Figure 29 with a focus on the root migration of the first 
mode. As can be seen, the first mode is stable throughout these dynamic pressures. However, a zoomed-in 
version of Figure 29 with a focus on the root migration of the third mode, presented in Figure 31, indicates 
that initially, the third mode exhibits a slight instability before becoming more stable. 
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Figure 30. Viscous FUN3D ROM root locus plot at M=1.141, zoomed-in to the first mode. 


24 of 29 


American Institute of Aeronautics and Astronautics 


320 - 


315 


310 


i_ 

vs 
c 

w 
vs 

E 305 


300 


- 0.2 


-0.15 


- 0.1 


-0.05 0 

Real, rad/sec 


0.05 


0.1 


0.15 


Figure 31. Viscous FUN3D ROM root locus plot at M=1.141, zoomed-in to the third mode. 


The four generalized coordinates from a FUN3D full viscous solution are presented as Figure 32. At 
this low dynamic pressure, all four generalized coordinates are lowly damped. Figure 33 and Figure 34 
present the first and third generalized coordinates from Figure 32, respectively. Visual analysis of both of 
these generalized coordinates indicates that these generalized coordinates appear to be stable although lowly 
damped. However, it is important to state a fundamental and important difference between a root locus plot 
and the visual, or otherwise post-processed, analysis of generalized coordinates over a short period of time. 
A root locus plot, by definition, exhibits the roots of a system as time approaches infinity or as the system 
reaches steady state. In contrast, the analysis of the initial transient response of a generalized coordinate 
over a short period of time can be deceiving as the response can change if the response was viewed (or 
analyzed) over a longer period of time. Therefore, based on the root locus plot at this condition, it appears 
that the third mode is unstable although that is not apparent in Figure 34. 
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Figure 32. Generalized coordinates from viscous FUN3D full solution at M=1.141 and Q=6 psf. 



Figure 33. First generalized coordinate from viscous FUN3D full solution at M=1.141 and Q=6 psf. 
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Figure 34. Third generalized coordinate from viscous FUN3D full solution at M=1.141 and Q=6 psf. 


V. Concluding Remarks 

A comparison of linear, inviscid, and viscous aeroelastic solutions for the AGARD 445.6 wing were 
presented. The linear solutions were generated using the linear equations available within the CAP-TSD 
code. The inviscid and viscous solutions were generated using the FUN3D code. The results presented 
consisted of two types of solutions: full CFD solutions whereby the CFD aeroelastic solution was computed 
in the standard time iterative approach and reduced-order modeling (ROM) solutions whereby unsteady 
aerodynamic and aeroelastic ROMs were generated using FUN3D and used to generate aeroelastic root locus 
plots as a function of dynamic pressure. An important conclusion from this research is that visualization 
of the aeroelastic root locus plots provided by the ROM approach enables a direct and more comprehensive 
interpretation of the aeroelastic behavior. Case in point, for the inviscid (Euler) solutions, the third mode is 
always unstable while the flutter instability associated with the first mode does not occur until about 300 psf. 
If this third mode instability is present in all inviscid solutions published by other researchers is not known 
at this point in time. However, the importance of being able to view a root locus plot in addition to time 
domain responses of the generalized coordinates was well established with this particular result. Another 
related conclusion to be made is that the standard approach of viewing (and analyzing) a short time history 
of the initial transient response of generalized coordinates may not be sufficiently accurate since it is possible 
that the steady state response has not been reached over this short time interval. Here again, the use of 
the root locus plots generated using the ROM approach enable the complete observation of the aeroelastic 
response since, by definition, the root locus plots represent the dynamics of the system at the final steady 
state of the system. 
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